#pragma once
#include "copy.h"
#include <iostream>

/**
 * Wrapper around cublasSsymv().
 * performs Z = alpha*A*X + beta*Y.
 *
 * @param uplo: whether or not this is upper or lower storage, 1 is upper, 0 is lower
 * @param alpha: scaling factor for the matrix A
 * @param A: symmetric matrix supplied in triangular format
 * @param X: vector to be multiplied by A
 * @param beta: scaling factor for Y
 * @param Y: vector to be added to A*X
 * @param Z: preallocated buffor for result
 */
void ssymv
(
    int _uplo,
    float alpha,
    Matrix A,
    Vector X,
    float beta,
    Vector Y,
    Vector Z
)
{
    int n, lda, incx, incy, size_y;
    char uplo;

    lda = (A.stride == 1) ? A.sub.stride : A.stride;
    n = A.length;
    uplo = _uplo ? 'u' : 'l' ;
    incx = X.stride;
    incy = Y.stride;
    
    if(n != X.length && n != Y.length){ //dimension mismatch
        std::cerr << "ssymv: ERROR mismatchin in the dimensions of X or Y and A" << std::endl;
    }

    size_y = (1 + (Y.length * incy)) * sizeof(float);
    cudaMemcpy(Z.data, Y.data + Y.offset, size_y, cudaMemcpyDeviceToDevice);

    cublasSsymv(uplo, n, alpha, A.sub.data + A.sub.offset, lda, X.data + X.offset, incx, beta, Z.data, incy);
}

void ssymv
(
    int _uplo,
    float alpha,
    Matrix A,
    stored_sequence<float> X,
    float beta,
    stored_sequence<float> Y,
    stored_sequence<float> Z
)
{
    Vector x(X.length, 1, X.data);
    Vector y(Y.length, 1, Y.data);
    Vector z(Z.length, 1, Z.data);
    ssymv(_uplo, alpha, A, x, beta, y, z);
}
